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ABSTRACT 

We present the results of a study of the low-mass X-ray binary 4U 1636-536. We have performed temporal analysis of all available 
RXTE/ASM, Swift/BAT and MAXI data. We have confirmed the previously discovered quasi-periodicity of =; 45 d present during 
~2004, however we found it continued to 2006. At other epochs, the quasi-periodicity is only transient, and the quasi-period, if present, 
drifts. We have then applied a time-dependent accretion disc model to the interval with the significant X-ray quasi-periodicity. For our 
best model, the period and the amplitude of the theoretical light curve agree well with that observed. The modelled quasi-periodicity 
is due to the hydrogen thermal-ionization instability occurring in outer regions of the accretion disc. The model parameters are the 
average mass accretion rate (estimated from the light curves), and the accretion disc viscosity parameters, a, for the hot and cold 
phases. Our best model gives relatively low values of a- co id ^ 0.01 and tthot — 0.03. 
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1. Introduction 


4U 1636-536 is a low-mass X-ray binary (LMXB) discovered 
by Willmore et al. \ \91A\ . The photometry of the optical coun¬ 
terpart (V801 Ara) shows a short orbital period of 3.79 h ( |Giles| 
et al. 2002). The binary system consists of a late-type, low- 
mass (- 0.3-0.4 Mg) donor, which transfers mass onto a neu¬ 
tron star (Fujimoto & Taam|1986(|van Paradijs et al.|1990)l. Gal¬ 


loway et al. (2006) estimated the distance to 4U 1636-536 to 
he D - 6.0 ± 0.5 kpc from Eddington limited X-ray bursts, as¬ 
suming the neutron star mass of 1.4M 0 and the stellar radius 
of 10 km. According to Casares et al. ( j2006 i, the mass func¬ 
tion and mass ratio of 4U 1636-536 are f(M) - 0.76 + 0.47 Mg 
and M 2 /MNS — 0.21-0.34, respectively, where Mns is the mass 
of the neutron star and M 2 is the mass of the donor. They also 
estimated the inclination as i - 36°-60°. The binary is a persis¬ 
tent X-ray source, although it shows significant flux variations 
on both long and short time scales. On time scales of hours, its 
flux varies by a factor of ~2-3 (|Hoffman et al.||1977 Ohashi| 


et al.|1982| Breedon et al.| 1 986] Hasinger & van der Klis| 1 989). 

The presence of kHz quasi-periodic oscillations, which are also 
visible in the system during X-ray bursts, shows that the neu¬ 
tron star has been spun-up through accretion ( Zhang et al.|1996 
S trohmayerjl 999). 


4U 1636-536 has been monitored daily in the 1.3-12.2 keV 
energy range by the All Sky Monitor (ASM) on-board of the 
Rossi X-ray Timing Explorer ( RXTE ) from 1996 until 2011. Dur¬ 
ing the first four years of RXTE /ASM observations (1996-2000) 
the source count rate was relatively stable at ~ 15 cts s 1 . After 
2000, it started to gradually decline and occasionally show a sta¬ 
tistically significant quasi-periodic variability ( Shih et al.|2005| >. 
Those authors reported the presence of a long-period, =* 47 d, 


Year 

1996 1998 2000 2002 2004 2006 2008 2010 2012 



MJD (days) 

Fig. 1 . The 50-d average RXTE/ASM light curve of 4U 1636-536 in the 
energy range of 1.3-12.2 keV from January 1996 to November 2011. 


quasi-periodic variability in the 2004 light curve. They suggested 
that the observed flux variability is caused by the variability of 
the accretion flow related to X-ray irradiation of the disc. 

In our paper, we study the variability of 4U 1636-536 taking 
into account the currently available data from three X-ray mon¬ 
itors, spanning almost 20 years (1996-2014). We interpret the 
data in terms of a disc instability model. In Sect. [2j we describe 
the X-ray data and perform their timing analysis. In Sect. [3] we 
present the theoretical model of evolution of an accretion disc 
around a neutron star used in the paper. We model the evolution 
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Fig. 2. The power spectra (with offsets) for the RXTE /ASM (left panel), the SwiftfBAT (middle panel) and the MAXI (right panel) in the time 
intervals given in MJD by: (A) 50088-50741, (B) 50742-51452, (C) 51453-52162, (D) 52163-52886, (E) 52887-53564, (F) 53565-54185, (G) 
54186-54749, (H) 54750-55322, (I) 55323-55869, (J) 55893-56461, (K) 56462-56738. 


of the accretion disc for a wide range of viscosity parameters. We 
find that the observed quasi-periodicity can be explained by the 
thermal instability due to the ionisation of hydrogen. We discuss 
our results, and finally, in Sect. [4] we give our main conclusions. 


2. Data analysis 


We analyse the X-ray data obtained with the RXTE/ ASM (Bradt 
|et al.|l 99~3]|Levine et al.|1996 l, the Burst Alert Telescope (BAT) 


on board of Swift (Mark wardt et al.|2005)|Krimm et al.|20l3] l and 

the Monitor of All-sky X-ray Image (MAXI) on board of Inter¬ 
national Space Station (Matsuoka et al. |2009| ). One-day average 
data from these instruments have been used to analyse the flux 
variability at different epochs. The ASM operated in the energy 
range of 1.3-12.2 keV with the sensitivity of ~ 30mCrab. The 
BAT observes 88% of the sky each day with the sensitivity of 
~ 5 mCrab in the energy range of 15-50keV. The SwiftfBAT has 
observed 4U 1636-536 since 2005. The Japanese X-ray camera 
MAXI has observed X-ray sources since August 2009, scanning 
most of the sky every 96 minutes. It operates in the energy range 
of 2-20 keV with the sensitivity of 20 mCrab. 

The X-ray flux history in the ASM data is shown in Fig. [T] 
The light curve can be divided into three main parts. From 1996 
to 2000, the flux was relatively stable (~ 15 cts s~'), with only 
some low-amplitude variability. Between 2000 and 2002, the 
flux gradually declined to a level below ~ 10 cts s 1 . Since 2003, 
the flux was below ~ 10 cts s and high-amplitude variability 
was present. 

We have linearly detrended all the data and then applied a 
Fourier transform. We use the standard Lomb-Scargle method 
(Scargle 1982), as implemented in the astroML (Vand erPlas] 
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Fig. 3. Top: one-day average RXTE /ASM light curve for 2004. Bottom: 
the corresponding power spectrum. 


|et al.|20f2) Python library. This method is suitable for our un¬ 
evenly sampled time series. We have divided all the data into 
11 intervals with the length of ^ 2 yr each, denoting them with 
the capital letters: A: 1996-1997, B: 1997-1999, C: 1999-2001, 
D: 2001-2003, E: 2003-2005, F: 2005-2007, G: 2007-2008, H: 
2008-2010,1: 2010-2011, J: 2011-2013, and K: 2013-2014. The 
power spectra of these intervals for the data obtained with three 
instruments are shown in Fig. [2] 

We then focus on 2-yr intervals from the ASM data shown 
in the Fig. [2] We see no periodicity in the interval A. From B 
to D, some quasi-periodicity appears, but is not statistically sig- 
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Fig. 4. The same as in Fig.[5]but for 2005. 
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Fig. 5. The same as in Fig.[3]but for 2006. 
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Fig. 6. The same as in Fig.[3]but for 2007. 
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Fig. 7. Top: one-day average ASM 1.3-12.2 keV (thin black curve) and 
BAT 15-50 keV (thick red line) light curve of 4U 1636-536 from MJD 
53755 to 53963 in Crab units. Bottom: sum of one-day average ASM 
and BAT light curves (1.3-50keV). 
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In order to perform a more detailed analysis, for the time 
span where significant quasi-periodicity is visible, we focus on 
the power spectra of the ASM data from the second half of the in¬ 
terval E to the first half of the interval G, i.e. 2004—2007. There¬ 
fore, for this time span the data were divided into 1-yr intervals. 
The light curves and the corresponding power spectra are shown 
in Figs. w Relatively strong quasi-periodic variability is vis¬ 
ible in 2004, 2005 and 2006, where the power spectra maxima 
are at P — 44.5 d, 44.0 d and 45.5 d, respectively. In 2007, the 
Lomb-Scargle power significantly decreased to P\s ~ 0.25 but 
also the power-spectrum maximum shifted to around 37 days, 
see Fig. [3] 

The BAT data are relatively noisy. We see some quasi¬ 
periodicity with P si 45 d in the interval F. Then, it seems that the 
period of the variability is decreasing until the interval I, where 
the 45-d periodicity is again visible. 


In order to approximate the bolometric light curve of 
4U 1636-536, the fluxes of the ASM and BAT in the Crab unit 
were summed. We show the individual light curves for soft and 
hard X-ray ranges and their sum in Fig. [7] for the time span of 
MJD 53755-53963 within the interval F. These are the most pre¬ 
dominant periodic data in our whole sample and therefore this 
part of the light curve was later used as a template for our the¬ 
oretical model. The apparent delay of the ASM light curve with 
respect to that of the BAT of about 10 d (calculated by cross¬ 
correlation) appears to be related to an anti-correlation between 
these bands reflecting the spectral softening from the hard spec¬ 
tral state to the soft one with the pivot energy of ~ 10 keV, see 
e.g. Done et al. ( j2007| l and references therein. In an ideal case, 
the maximum of the flux above the pivot would be at the mini¬ 
mum flux below the pivot, which corresponds to the delay of a 
half of the period. In our case the pivot is within the ASM energy 
range, which then reduces the delay. 


nificant, with the peaks in the power spectra being low. Finally, 
from the interval E (2003-2005), significant quasi-periodicity is 
visible. A similar periodicity was noticed by Shih et al. ( |2005[ > 
for the 2004 epoch. The quasi-periodic variability began to grad¬ 
ually decline from the interval G. 


We see some ~40-d quasi-periodicity in the first 2-yr inter¬ 
val, H, of the MAXI data. In the next interval, I, there is a strong 
quasi-periodicity with P ^ 45 d, similarly to the corresponding 
Swift/HAT data. Later, this periodicity disappears and another 
one with approximately doubled period, i.e. ^70 d, appears. 
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3. The theoretical model 

The theory of geometrically thin, stationary, accretion disc of 


Shakura & Sunyaev (19731 describes well accretion as the 


source of energy of many astrophysical objects. On the other 
hand, many of the observed accretion sources are strongly vari¬ 
able, which may be related to non-stationary effects predicted 
theoretically. In particular, [Meyer & Meyer-Hofmeister| ( |1982] > 
and |Smak| ( |1982) > independently found that accretion discs be¬ 
come unstable in the range of temperatures corresponding to ion¬ 
ization of hydrogen. This instability is the result of an inverse re¬ 
lation between the temperature and opacity in the range between 
~ 10 3 K and ~ 10 4 K. The Rosseland mean opacity is smoothly 
decreasing with increasing temperature at higher temperatures. 
But when the temperature drops below 10 4 K, hydrogen partially 
recombines, and the opacity steeply decreases with decrease of 
temperature. Then, at T < 10 3 K, hydrogen becomes completely 
neutral, and the huge grain opacity dominates ( |Alexander|1975[ >. 

The theory of partial hydrogen ionisation instability was 
originally developed to explain the large-amplitude luminos¬ 


ity variations of cataclysmic variables (dwarf novae; Meyer & 
|Meyer-Hofmeister|1982[|Smak|1984| ). It is also believed that the 


same mechanism is responsible for outbursts in soft X-ray tran- 

sients 

Cannizzo et al.||1982| Dubus et aL||2001ULasota 

2001 f 

Bagiris 

ca et al.|2014) and active galactic nuclei (|Lin & S 

fields 

1986 

Clarke 1988; Mineshige & Shields 1990; Siemiginowska 

et al. 

1996). 


Also, Hameury et al. ( 1998) 1 studied time dependent numer- 
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We use the method of |Smak[ ( |1984| i in order to calculate the 
time evolution of an accretion disc around a neutron star. We 
follow the changes of the disc temperature and density. Then the 
surface temperature determines the emitted luminosity, observ¬ 
able in the X-ray band. We do not explicitly take into account 
the emission of the boundary layer, but include it in the estimate 
of the total luminosity as a function of the rate of mass supply at 
the disc outer radius. 

3. 1. Vertical structure of an accretion disc 

In a geometrically thin accretion disc, the vertical and radial 
structures are decoupled ( [Pringlejl98l| . The equations that de¬ 
scribe the disc vertical structure are very similar to those describ¬ 
ing the internal structure of a star. For a Keplerian disc with the 
angular velocity, Q = ( GM/R } j 1 ^ 2 , we have the energy balance. 


d F dQ 

— = -aP — , 
d- d R 


(1) 


where M is the mass of the accretor, F is the energy flux, a is 
the viscosity parameter, P is the pressure, and z and R are the 
height and radius of the accretion disc, respectively. Then, the 
hydrostatic equilibrium reads. 


1 dP 

--r- = ~^z = g z , 

p dz 


( 2 ) 


where g z is the vertical component of gravity. The energy trans¬ 
fer equation can be written as 


ical models of the accretion disc evolution to explain the cyclic 
outbursts of dwarf novae and soft X-ray transients. They used 
an adaptive grid technique, in which the outer radius of the disc, 
P 0 ut, is not fixed but it varies with time. In their simulations, the 
variable R out allows for a propagation of the heat front from outer 
disc radii inwards, which is aimed to explain the properties of the 
dwarf nova SS Cygni. 

Similarly, |Dubus et al. ( 2001| > used the adaptive grid numer- 


d In T dlnP 


d; 


d~ 


(3) 


where T, P and p are temperature, pressure and density, respec¬ 
tively, and V s dlnT/dlnP. In an optically-thick disc without 
convection, it is the radiative gradient, 


V - V la d - 


ical scheme. A large number of grid points, at least about 100- 
800, is required in such computations in order to conserve the 
mass of the disc, which was checked a posteriori by integration 
of the surface density over radius. In our computations, the outer 
disc radius is fixed, and the mass is conserved. We use a linear 
grid in 2 \fr. We keep the number of grid points at 150, to allow 
for moderate resolution. In this way, we assure that the thick¬ 
ness of the disc is not larger than the size of the disc cell, A r, 
so the assumption of the geometrically thin disc is satisfied. This 
condition is important because we calculate the 1 + 1D structure 
(see below) rather than perform 2D hydrodynamical simulations. 
Still, when we increase the number of grid points to, e.g., 300, 
the obtained light curve is only moderately affected, keeping the 
qualitative character of the variability. 

[Lipunova & Shakura| ( |2000[ |2008| > presented an analytical 
model of time dependent accretion disc, in which the continuity 
equation for the surface density change was solved. In this way, 
the authors explained the power-law decay of the luminosity and 
model the light curves of the X-ray novae. These authors do not 
compute numerically the thermal-viscous evolution and cyclic 
outbursts of the sources, because this would require solving the 
time-dependent thermal balance. Nevertheless, Sulei manovetal.| 
(2008) successfully explained with such a model the decay pro- 


kPF 

4FradCg, 


(4) 


where k is the Rosseland mean opacity and P rM [ is the radiation 
pressure. However, in the case of partially ionized hydrogen, a 
major mode of energy transport may be due to convection, V conv , 
which acts together with the radiative transport. We use here V = 
V CO nv + V la d, where V conv is calculated using the mixing-length 
theory, as in stellar models ( Paczyriski|l969| l. 

The equations are integrated between the disc midplane, 
Z — 0, and the disc photosphere at the height II, given by the 

poo 

optical depth of the top layer of L pxdz = 2/3, where p is the 
mass density (see |Rozariska et al.||1999). The tables of opac¬ 
ities are taken from |Cox & Stewart| ( |1969| |1970| ) and |Seaton| 


et al. (1994l for temperatures higher than log 10 T > 3.5 and 


from jAlexander et al.| l fl98'3) for lower temperatures. At the pho¬ 
tosphere, the effective temperature of the disc is T e g = T(H), 
and <x7/ 4 (| = Q - Q is resulting from the heating of the disc 
through the viscous energy dissipation, Q + , balanced by the ra¬ 
diative cooling, Q~, where <x is the Stefan-Boltzmann constant. 
The integration of equations for the vertical structure gives the 
surface density, E, 


E = 2 


/' 


pdz. 


(5) 


files of the black-hole transients A 0620-00 and GX 1124-68. In 
our computations, we solve both the time-dependent evolution of 
the disc density and temperature; thus, modeling of cyclic out¬ 
bursts and quiescent phases is possible. 


The viscous energy dissipation rate per unit time and volume 


Q+ = -aP to tC2, 


( 6 ) 
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Fig. 8. An example of the surface density-effective temperature E-Jeff 
relations (S-curves) for an accretion disc around a 1.4 M G neutron star 
at R = 2 x 10 4 R g (solid black and red dashed curves) and 2 x 10 5 R g 
(thick blue curve). The black curve corresponds to a co i d = 0.01, ar hot = 
0.05, the red dashed curve, to ar co i d = 0.01, a^ot = 0.03, and the thick 
blue curve, to a cold = a hot = 0.01. For the black curve, the hydrogen 
ionisation instability begins in A and ends at B: analogous points can 
be determined for other curves. 


where we assumed that the viscous stress tensor R-<p component 
is proportional to the total pressure P lM (sum of the gas and radi¬ 
ation pressures), and scaling with a ( [Shakura & Sunyaev|1973| l. 
Since we know the relation between the effective and central 
temperatures from the energy transfer equation, we can repre¬ 
sent the energy balance in a X-7i. (l ' diagram. The 7’ c n is directly 
related to the accretion rate, M, which average is given by exter¬ 
nal conditions (i.e., the mass transfer rate from the companion 
star) and which is a parameter of the model. Example solutions 
are plotted in Fig. [8]for the accretor mass of Mns = 1.4 M Q and 
selected values of the viscosity parameter, a, and the radius, R. 

We can distinguish three characteristic regions in the S- 
curve. The lower branch corresponds to thermally stable and 
cool neutral hydrogen. The upper branch corresponds to ther¬ 
mally stable but hot gas, where hydrogen is fully ionised. The 
opacity is then due to bound-free transitions in heavy elements, 
free-free transitions and electron scattering. Finally, the middle 
branch corresponds to partially-ionised, thermally unstable re¬ 
gion of the disc ( Bath & Pringle|1982||Faulkner et al.|1983| >. 

Calculating the time-dependent radial structure of an accre¬ 
tion disc can be significantly simplified because of the universal¬ 
ity of the above properties. The critical column densities, X mm 
and Z max , between which the unstable branch of the stability 
curve appears (i.e. the points A and B in Fig. [8]), are functions 
of the accretor mass (fixed in our case), the viscosity parame¬ 
ter, and the radius. On the other hand, the corresponding critical 
values of the effective temperature are given by the temperatures 
at which hydrogen becomes dominantly neutral or ionized, and 
thus are nearly universal, i.e., only weakly depend on the accre¬ 
tion flow parameters. 

Smak| (1984) suggested the following simple scaling rela¬ 


tions, assuming that the shape of the stability curve does not 
change and the critical points shift with the radius throughout 
the disc. 


log ^ = log J e B ff + a, log(R/2 x 10 4 /( g ) + a a log(a/0.1), (7) 

log X A = log E b + b r log(R/2 x 10 4 /C) + b a log(or/0.1), (8) 


where 2 x 10 4 R g is a chosen reference radius and R g = GM/c 2 is 
the gravitational radius. In our procedure, we choose two pairs of 
values of R and a , for which we calculate the Z-7 c fr curves. We 
then calculate the positions of the A and B points, and then fit 
the coefficients to the stability curves calculated at intermediate 
values. The resulting values are a r = -0.09, a„ = -0.01, b r = 
1.11, b a = -0.73. 

In order to reproduce the observed amplitudes of dwarf no¬ 
vae outbursts, it was suggested that the viscosity parameter a 
must have different values on the upper and lower branches 
of stability curve (see, e.g., Lasota |2001 Done et al. 2007 1 
for reviews). In our computations, we also assume it, implying 
£min = £ A - S(ffhot) and E max = I B = 2(a cold ). A bridging for¬ 
mula is used for the transition between the two values of a along 
the stability curve. 


l°gio( a ) = logioOcoid) + 


^ 4 g I O^hot Arcoid) 

1 +(2.5 x 10 4 K/r c ) s 


(9) 


The above old, observationally motivated, approach, was 
confirmed by numerical simulations of the magneto-rotational 
instability (MRI, |Balbus & Hawley11991) . This instability is sup¬ 
posedly the physical process that is acting behind the viscosity 
mechanism in astrophysical flows. The non-linear development 
of the MRI leads to the disc turbulence, and hence turbulent vis¬ 
cosity, because magnetic fields are frozen into the fluid in a dif¬ 
ferentially rotating disc. Thus, the MRI provides an efficient way 
to transport the angular momentum outwards in the disc and al¬ 
low for accretion. It has been shown, e.g., in laniuk et al. ( |2Q04] >, 
that the resistive diffusion of magnetic field, as quantified by the 
magnetic Reynolds number, is low in a cold state of the discs 
of stellar mass accreting compact objects (white dwarfs, neutron 
stars and black holes in binaries). Therefore, it is expected that 
the MRI turbulence disappears in the quiescent disc. In the out¬ 
burst state, on the other hand, the action of MRI turbulence is ef¬ 
ficient. Consequently, the approach postulated for these systems, 
namely the cr co id < ahot is justified, and we follow it here. 

In our numerical approach, we first calculate the disc verti¬ 
cal structure as a sequence of local solutions. We then calculate 
the time-dependent radial structure, which yields the final light 
curve. This method is described in detail by |Smak ( 1 984[ > and 
Siemiginowska~et al.|(|1996|). 


3.2. Results and discussion 


Following Galloway et al. ( |2006| ), we assume Mns = 1.4Mo, 
7?ns = 10 km, and D = 6.0kpc. We have estimated the aver¬ 
age bolometric flux by scaling the ASM and BAT spectra to the 
Crab, adding them and then applying a correction for the emis¬ 
sion beyond the energy ranges of those instruments of 50%. This 
yields Fboi - 4.4 x 10~ 9 erg cm -2 s _1 , and, assuming isotropy, 
Lboi = 4nD 2 F ho i - 1.9 x 10 37 erg s 1 . The accretion efficiency in 
the Newtonian approximation (e.g., Frank et al.|2002)> is 


1 = 


GMns 

Rnsc 2 


= 0 . 21 , 


GO) 


which implies the accretion rate averaged over an outburst of 


M=^~1.6x 10~ 9 M 0 yr _ 
rj c z 


(id 


We assume M — 1.5 x 10 9 M G yr in the calculations. We study 
models with 0.001 < a co id <0.1 and with the ratio of 0.14 < 
a+oid/ffhot < 1.0. We assume the disc outer radius equals R out = 
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Fig. 9. The model 1 for or cold = 0.01, a hol = 0.05 (the two upper panels), 
and model 2 for a cold = 0.01, a hot = 0.03 (the two lower panels). For 
each model, we show the light curve and the corresponding power spec¬ 
trum in the upper and lower panel, respectively. The solid curves show 
the observed light curve, given by the sum of one-day average ASM 
and BAT fluxes, with the time of zero corresponding to MJD 53755. 
The dashed curves show the theoretical light curves and their power 
spectra. 


2 x l0 5 R g , which corresponds to about 80 per cent of the Roche- 
lobe radius of the neutron star. The disc inner radius is at 6 R g . 

We have selected time span from 20 January 2006 to 16 
August 2006 for the ASM and the BAT light curves, and we 
summed the fluxes from them as shown in Fig. [7] This light 
curve shows relatively regular outbursts with high amplitude. As 
we found in Section[2] the quasi-period of this part of the data is 
P =* 45 d. We have calculated a large numbers of models, com¬ 
paring them to the data. 

We show here two models, 1 and 2, in Fig. [9] For each model, 
we show the light curve and the corresponding power spectrum 
in the upper and lower panel, respectively. We see that our model 
1, with a co ]d = 0.01, ayoi = 0.05, reproduces correctly the ob¬ 
served outburst period, but it significantly overestimates its am¬ 
plitude. Then, our model 2, with a- co | d = 0.01 and ah 0 t = 0.03, 
approximately reproduces both the period and its amplitude. In 
this model, the instability begins at > 0.5 R out . 

We find that the outburst amplitude decreases with the in¬ 
creasing acoid/aTot- Then, a^oid determines the interval between 
outbursts and cthot, the time scale of a single outburst. Both time 
scales decrease with the increasing value of a. We assume here 
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is the constancy of the average M, which is justified by the con¬ 
stancy of the flux averaged over an outburst in the chosen data 
interval. However, changes of the M averaged on time scales 
longer than that of the outburst will lead to changes of the quasi¬ 
period and the amplitude of an individual outburst, as seen in the 
data. Specifically, for our chosen parameters, an increase of M 
above 1.5 x 10 9 M Q yr 1 reduces the outburst amplitude and the 
duration of the minimum, while it increases the duration of the 
maximum. A decrease of M leads to the opposite behaviour. The 
instability disappears at M 5 x 10~ 9 M G yr 1 . Thus, our model 
can explain the overall properties of the entire data set. 

We note that Shih et al. ( 20051 proposed that the observed 
quasi-periodicity in 4U 1636-536 is explained by an instabil¬ 
ity due to the irradiation of the disc by X-rays. However, they 
did not perform any calculations to support their proposal. Our 
model does not invoke irradiation. We note that irradiation, due 
to the additional heating imposed on the disc by the X-rays, has 
a stabilizing effect. Therefore, it may still change the quantita¬ 
tive details of the disc instability model ( |Lasota|2015| ). As it was 
shown, e.g., by Dubus et al. ( 2001| >, somewhat longer decay- 
phase and recurrence Times are obtained if irradiation is taken 
into account and the inner disc is allowed to be truncated be¬ 
tween the outbursts. In our calculations, we do not take the disc 
truncation and irradiation into account. 

The different values of a in the cold and hot states that we 
found here to be necessary, are in general agreement with these 
suggested by observations of dwarf novae and soft X-ray tran¬ 
sients. As it was found by |Hirose et al.| ( j2014[ ), an intermittent 
thermal convection on the upper stable branch of the disc sta¬ 
bility curve may strengthen the magnetic turbulence and sig¬ 
nificantly increase the value of a. Nevertheless, the values of 
a found here are relatively low, in the range of 0.01-0.03 for 
our best model, and the increase of its value on the hot branch 
is relatively modest. This means that the values of the Maxwell 
and Reynolds stresses in the cold and hot disc states are rela¬ 
tively similar, as it was also found in some shearing-box simula¬ 
tions (see [King et al.|2O07| and references therein). Thus, at least 
for the neutron star binary studied here, the poloidal magnetic 
field required to enhance the angular momentum transport in the 
ionised disc cannot be very large, and the local MRI turbulence 
is sufficient to account for its behaviour. In particular, the results 
obtained here imply the value of /3, the gas to magnetic pressure 
ratio, to be on the order of ~ 1/(2 a) (Yuan & Narayan[ |2014| l, 
which corresponds to ft ~ 15-50. This gives tighter constraints 
than those resulting from other numerical simulations (see | Yuan] 
& Narayan 2014 for a review). 


4. Conclusions 


We have analysed all the currently available X-ray data of 
LMXB 4U 1636-536 from the RXTE/ ASM, Swift/BAT and 
MAXI detectors, spanning together the energy range of 1.3- 
50 keV, and observing the source from 1996 to 2014. We have 
confirmed the previous finding of an appearance of the quasi¬ 
periodicity in the X-rays from the source in 2004 ( Shih et akj 
|2005| and found it continued until 2006. The average period 
during that epoch is -45 d. However, the quasi-period drifts (in¬ 
creasing up to ~70 d or decreasing up to -37 d) or disappears 
altogether, and its amplitude changes after 2006. The period of 
-45-d appears again during 2010-2011. 

We then have considered a model of the hydrogen-ionisation 
disc instability. Using the formalism proposed by |Smak ( |1984} > 
to explain outbursts in cataclysmic variables, we have found that 
it can be successfully used for this accreting neutron-star. In or- 











































Mateusz Wisniewicz et al.: Long-term quasi-periodicity of 4U 1636-536 resulting from accretion disc instability 


der to reproduce a part of the X-ray light curve of 4U 1636-536 
with a stable quasi-periodicity, we have numerically determined 
the vertical structure and calculated the radial time-dependent 
evolution of the accretion disc. We have estimated the average 
accretion rate based on the source distance, the X-ray light curve, 
and the accretion efficiency for the assumed mass and radius of 
the neutron star. The main free parameters of our model are then 
the values of the viscosity parameter, a. 

We have found that, in order to reproduce the observed light 
curve given the source parameters, we need to use two differ¬ 
ent values of a for the cold (with mostly neutral hydrogen) 
and hot (with ionized hydrogen) phases, which confirms a num¬ 
ber of previous findings, see, e.g., Lasota (2001]), [Done et~aLl 
(2007). Among the large number of the calculated models, the 
model with relatively low viscosity parameters, a co id = 0.01, 
Q'hot = 0.03 best reproduces both the amplitude and the period of 
the used part of the X-ray light curve. 
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